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Introduction 


There are a few things needed to be mentioned about this Technical Report #2. Firstly, the last 
section of this report contains what was to be included in what would have been Technical Report 
#1 Version 2. This means that there will no longer be a new version for Technical Report #1, 
and for other technical reports as well. Instead, the next technical reports will built upon their 
predecessors, as well as upon other things. 

By doing it this way I hope that the organization of the reports will become even simpler. The idea 
is to link research works and studies together so that they may be easier to find when needed. 
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Second order singularly perturbed sytem 


2"4 May 1998 
Consider the system described by a differential equation 


da: da 
e—~ + — =4, 1 
de®’ at ) 
or written as state equation 
1 = 9 (2 
Elo = —Xo+u. (3 


Some knowledge about the system 


By letting the right-hand side of the state equation be zero, every point on the 2,-axis is an 
equilibrium point. 
Equation 2 and 3 represent a model in standard form |Kha96] since Equation 3 has a real isolated 
root when ¢ = 0. 
(A model 
bf Gee) ee RE (4) 
Bo. .= Ota, ze). oe Re” (5) 


is said to be in standard form if and only if [Kha96] 
= g(t, Z, Z, 0) 


has k isolated real roots 
22h e),. Has ak 


And the model in Equation 4 and 5 reduces to a quasi-steady-state model or a slow model 


ef a the) 0 se (6) 


= 1 
x=(2x1,0) 0 — é 


the Eigenvalues are 0, —2. Thus the system has an equilibrium subspace and the qualitative be- 
haviour of the trajectories is summarized in Table 1. 
From Table 1, the system is stable when ¢ > 0 and unstable when ¢ < 0. 


) 


From the Jacobian matrix 
OF 


Aa ae 
Ox 


Response with various values of ¢€ 


Time responses when u = 0 of the state variables of this system were shown in Figure 1 (a) to (f), 
and the state plane were shown in Figure 2 (a) to (f). Figure 1 and 2 shows that the system is 
stable when ¢ > 0 and unstable when ¢ < 0. When |e] becomes smaller the response of the system 
becomes more rapid as 


Lo = -V 
E 


Technical report #2, presented to Professor kK. Furuta by Kittisak Tiyapan 3 


Values of € Nature of trajectories 
e>0 All trajectories converge to the equilibrium subspace 
e== 0 Degenerates, reduced to a 1*‘-order system 
Ee <0 All trajectories diverge away from the equilibrium subspace 


Table 1: Qualitative behaviour of the trajectory. 


rapidly coverges to its root x2 = 0 when ¢ > 0. 
The closer ¢ gets to zero the more parallel the trajectory becomes to the x2-axis. Equation 1 when 
€ + 0 and u = 0 would become 

z,=0, or r=0, 


which has the solution 
xi(t) =2,(0), and x(t) =0. 


The singular point ¢ = 0 is the point of discontinuity because the solution is stable when ¢ > 0* 
and unstable when ¢ > 0-. 
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(e) (f) 


Figure 1: Time responses of state variables. (a) ¢ = —1, (b) « = —0.1, (c) e = —0.03, (d) « = 0.01, 
(6/605 and (7) 2 — 1. 
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Figure 2: State planes. (a) ¢ = —1, (b) « = —0.1, (c) e = —0.03, (d) ¢ = 0.01, (e) ¢ = 0.1 and (f) 


aa 


Furuta Lab., Tokyo Institute of Technology 


The phase planes for various € are shown again in Figure 3, this time with the Initial Condition Set 


1 as the initial conditions. When ¢ gets smaller the slope of the trajectory becomes steeper as the 
rate in which x2 approaches its root becomes faster. 


Initial Condition Set 1 The set of n points described by 


n=0 

TOK 1 = $5... HS gob 
#1(O): S° a5, xe O(0) 4 
x1(0) = i, x2(0) = -i 
n = nt+2 

endfor 


29 order singularly perturbed system 24 order singularly perturbed system 
ularly ps ys ularly p 5 
=5 e=1 


(a) 
2" order singularly perturbed system 2" order singularly perturbed system 
e=0.1 £= 0.001 
5 ¢ . “ 5 
| \ 
at \ \ at 
\ 
3 \ © ? 3 
| | | 
2b \ \ | 2b 
tL | | 
\ | \ 
1 | \ ) ry \ | 1 
| | | | 
Sok x x x x ' j x * x “0 x 
| \ \ \ | 
fe 4 \ ® b \ \ 1 
| \ | 
\ | | | 
-2r | \ \ \ -2b 
\ \ \ \ 
\ \ | 
3 \ b ry \ 3b 
| \ 
\ | 
| \ ab 
| | 
6 4 2 0 2 4 6 6 4 2 0 2 4 6 
x x 


Figure 3: State planes (a) ¢ =5, (b) e=1, (c)e =0.1, (d) ¢ = 0.001 
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With feedback of a sign function 
Feed the sign function of the states 
u = sgn(x, + X2) 
back into the input and the results was obtained in Figure 4 (a) to (f) and in Figure 5 (a) to (f). 
(Though this may not be a good choice of input comparing with, for example u = K sgn (Cx; + 22), 
when C,K € R.) Figure 4 shows that the response of the controlled system changes more rapidly 
when |e| becomes smaller. The system is stable when ¢ > 0 and unstable when ¢ < 0. Both Figure 
1 in the previous section and 4 (d) to (f) show the boundary-layer characteristic of the system. 
In both figures the boundary-layer is t + |0,0.05], [0,0.4], and [0,1] when ¢ = 0.01, 0.1, and 1 
respectively. The increase in speed of the response as € decrease in magnitude applies to the case 
when € < 0 as well as when ¢ > 0. In other words, the response becomes faster as |e] — 0. 
The phase planes when the system is unstable (Figure 5 (a) to (c)) show only a small part of the 
trajectory due to the fact that the trajectory goes very far away from the origin. In all intent and 
purpose they only show that the system is unstable and the trajectory runs away. 
With the control as shown in Equation the Equation 2 and 3 become 
Ly = (7) 

efg = —£+sgn (41+ 22). (8) 

The root of Equation 8 when ¢ = 0 


0 = —x2 + sgn (41 + £2) 


Eo = sgn (Z, + Z2) 
which can not be written in a form 
To = lt): (9) 
In other words, Equation 9 has no isolated root. Therefore one can not isolate the fast mode from 
the slow mode by introducing a new variable 


Y =X%.—h(t,x7). 
But for the purpose of finding the boundary-layer model, let 
h(t, £1, %2) = sgn (a1 + 22). (10) 
Substitute x2 in Equation 8 with x2 obtained from 


Y = &—h(t,x1, 22) 


rq — sgn (41 + 22) 
tq = y+sgn (x1 +22). 


And introduce a new time variable tT obtained by 


dy dy 
ay _ ay 11 
“at ae ae 
dt 1 
(ns 12 
dt € (12) 
t 
T = TH] tdt 
to 
— ea 7% = 0 (13) 
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The new time variable 7 is generally known as a stretched time variable. 


Then 
d d 
€ a -- ee (v1 +2%2)|} = -y (14) 
d d 
oy -++ —sgn (xy + £2) = —Y (15) 
dr dt 
dys. a 
—— + —-sen (a1 + y+sgn (a1 +22)) = —y (16) 
dr at 
dy 


d 
ae —y — 4 sen (x1 +y+sen (21 + £2)) (17) 


Therefore the boundary-layer model can be represented by either one of the following, 
—y — £sgn (ex: ty +1), v1 +22 >0 


—y- “sen (x1 +y), 4, +%,2=0. (18) 
—yo ben (eit y— 1). wi Hey <0 


dy _ 
dr 


To find the reduced model look firstly at the three possible values which x72 may take 


+1, T1 +X. > 0 
tg= 0, %t+%=O0. (19) 
—-1l *%,4+%<0 


Then the reduced problem becomes 
+1 
=< 0 (20) 


with these three roots respectively. 
The discontinuous nature of the sign function makes the analysis of the variable structure control 
system difficult. This is due to the fact that one can not analytically find 


Osgna(t)  Osgn x(t) 
Ot é Ox; : 


where x(t) is an n-tuple vector. 
One way to overcome or go around this problem is to use a numerical approximation for the sign 
function where necessary. 
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Figure 4: Time responses of state variables. (a) ¢ = —1, (b) « = —0.1, (c) e = —0.03, (d) « = 0.01, 
(e) € =0.1 and (f) e =1. 
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Phase plane 
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Phase plane 


Figure 5: State planes. (a) ¢ = —1, (b) « = —0.1, (c) e = —0.03, (d) ¢ = 0.01, (e) ¢ = 0.1 and (f) 


e= 1. 
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With feedbacks of a sign function and state variables 
Feed the states as well as the sign function of the states back into the input 
u = —(2, + £2) + sgn (x1 + 22) (21) 


and the results was obtained in Figure 6 (a) to (f) (state variables) and Figure 7 (a) to (f) (phase 
planes). Similar to the previous studies, the system is unstable for ¢« < 0 and stable for ¢ > 0. 
Figure 6 (d) to (f) (€ > 0) shows that the system has a two-timescale property which is the 
characteristic of a singularly perturbed system. Having a fast response within the boundary-layer 
and a slow response elsewhere, the boundary-layer becomes narrower and the response faster as 
eé — 0*. Figure 6 (d) to (f)shows that the boundary-layer is t  [0,0.05], [0,0.23], and [0, 0.54] 
when € = 0.01, 0.1, and 1 respectively. 

The control input of Equation 21 turns the system (Equation 2 and 3) into 


Ly = 9 (22) 
éfg = —2, — 2%. +sen (x1 + 2X2). (23) 
Equation 23 when ¢ = 0 becomes 
1 1 
0 = —2, — 2%. +sgn(%, +22), or = al + 58gn (x, + £2). (24) 
In order to find the boundary-layer model, change the variable x2 to 
1 1 
Y= too) — to 521 — 5880 (x1 + 22). (25) 
in other words 
1 1 
2 = y-—~21+ =sgn (x1 + 22) (26) 
2 
1 1 
= y- 5f1 + 5sen (a, +y — =x, + =sgn (214+ %2)), (27) 


which is essentially 


i 1 1 1 
U— sa Seen (ye toi +4), Lite U 


T=) Yr al = 58g Yor }r1) ’ Ti+t2=0 . (28) 


1 1 1 1 
ym srit seen (yor hor — 3), Ci re CU 


And use the stretched time variable 7 in Equation 12. 
Substitute Equation 28 into Equation 23 to get the boundary-layer model 


dy ldzx 1 
E a. 5 3 588n (x1 +r 2) = —2y 
d 1 1 1 
€ i - ; | re _ zoe (a, +22) + 5oen (ay + 2) = —2y 
d 3 1 1 
. Spy = Gt ee (a, + 22), 
which is essentially 
—3y—19¢,-41, 2,+2.>0 
dy rd 4 17 gq 1 2 
ae ye ae a Li pas = 0 (29) 
Sy ae ae wie aS 0 
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Figure 6: Time responses of state variables. (a) ¢ = —1, (b) « = —0.1, (c) e = —0.08, 
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Time in seconds 


L 
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Time in seconds 


(f) 


(d) « =0.01, 


(e)€=0.1 and (f) e =1. With a sign function feedback and state variable feedbacks. 
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Figure 7: State planes. (a) « = —1, (b) e = —0.1, (c) « = —0.03, (d) « = 0.01, (e) ce = and (f) 
e=1. With a sign function feedback and state variable feedbacks. 
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Variable structure control 


Let the hyperplane be described by 


8 =cCrt,+ 2%, (30) 
which has the derivative 
a CL in Lo. (31) 
And let the input be 
u=—Ksens = —Ksen(cr; +22), K ERT (32) 
Choose a Lyapunov function as 
1 
= 53 (33) 
Recall that the state equations are 
Ly = (34) 
i K 
tg = ——@ — —sgn (cx, +22). (35) 
E E 
The derivative of the Lyapunov is then 
V = s6 (36) 
1 K 
= (cr, +22) (cxs — 582 — sen (cx, + 2) (37) 
3 i K 
= Car —- See us cx. — v2 en con (cr, + 22) — sen (cx, + £2) (38) 
E E 


To observe the effect of ¢ on the stability of the controlled system, let 
C=27,, “and: Jo 233, 


and plot V against x, and zz as shown in Figure 8. From Figure 8 whether V>0 or V < 0 depends 
on the value of ¢. For example, 4x such that V > 0 when « is 1.7 and 0.08, but V < 0 Vx when 
é = 0.3. 


Technical report #2, presented to Professor K. Furuta by Kittisak Tiyapan 15 


Figure 8: Derivative of Lyapunov, ¢ = 1.7 for (a) and (b), « = 0.3 for (c) and (d), and ¢ = 0.08 
for (e) and (f) 
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Figure 9 (a) and (b) shows the plots of the case of Figure 8 (c)-(d) and (e)-(f) respectively. Initial 
conditions were the Initial Conditions Set 1. Similar to Figure 9 (b), Figure 10 shows the case 
where ¢ = 0.08 but with the Initial Conditions Set 2 as the initial conditions instead. Despite the 
plot of V which shows that the system is unstable for ¢ = 0.08 the state planes of Figure 9 (b) and 
Figure 10 are stable. 


Initial Condition Set 2 The set of n points described by 


n = 0, where n is the number of initial conditions 
i= 8 
for j = -8, -5, -2, 1, 4, 7 

e100): S21 2200) S94 wn Set 

e100} =: =2y e200). = jf, ee 


B10) = 3),-3200) = 1, aS net 
e100). S35; 3200) == 1,, a = Fl 
endfor 
The original system (Equation 2 and 3) has eigenvalues at \; = 0 and Ay = —12.5, with the 
eigenvectors at 
2 sie —0.0797 
1=to }? | 0.9968 
for Ay and A»2 respectively. 
The system is linear since 
AK a, (39) 
ie A is a constant matrix. Also A is not a stability matrix since the condition 
Red; <0 


is not satisfied. The equilibrium point at the origin is not unique as it has already been pointed 
out earlier that every point on the x-axis is an equilibrium point. The x,-axis is a nontrivial null 
space of the matrix A 

Change of variables 


1 <=0,0797 | (40) 


0 0.9968 


leads to 


1 (9 0 Ge 0 1 
a Wns ease ali ee one Nei BS 
The solution to this system is 


31(t) = 21(0),  za(t) = z2(0)e7*?* = 20(0)e*. 
The fact that A is not a stability matrix can be seen by trying to solve the Lyapunov equation 
PA+A'TP==Q, (41) 


where Q is a positive definite symmetric matrix. The Lyapunov equation yields no solution for this 
system. @ 
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2" order singularly perturbed system 2" order singularly perturbed system 
€ = 0.3, with variable structure feedback € = 0.08, with variable structure feedback 


—______—_9 


Figure 9: State planes (a) ¢ = 0.3, (b) ¢ = 0.08 
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order singularly perturbed system 
€ = 0.08, with variable structure feedback 


-8 


-2 


Figure 10: State planes ¢ = 0.08 


Technical report #2, presented to Professor K. Furuta by Kittisak Tiyapan 19 


Figure 11 (a) to (f) shows some more of the plot of V with various parameters. Figure 11 (a) and 
(b) have 
C= 2.45. cand kh =38 


as in Figure 8 but here 
em 0. 


Figure 11 (c) to (f) have 

Kk = 100 
while ¢ = 1.7 and ¢ = 0.08 for (c)-(d) and (e)-(f) respectively. Figure 11 shows that when ¢ < 0 
the controlled system is unstable (Figure 11 (a)-(b)). For « > 0 a larger value of K gives a wider 
range of €¢ which the controlled system is stable. In other words the larger the value of gain K of 


the control input is, the more robust to € this control scheme will be. Compare Figure 8 (a), (b), 
(e), and (f) with Figure 11 (c), (d), (e), and (f) respectively. 


N.B. One way to stabilize the origin of the system is perhaps to solve the equation 
A’PE+ E"PA-(E™PB+8)R-'(B'PE+S")+Q=0, (42) 
or equivalently 
FT PH LE’ PFE PBR "BPE +Q—SR9' S0, (43) 


where F = A~BR7!S", Q and R are symmetric positive definite matrice, and the system is written 
in the form 
Ez = Ar+ Bu. 


As an example, let 
Oa hake = 


12.62 1 
pal 1 satel 


will lead to 


The feedback using this method is then 
u=—-Gs, G=R'(BTPE+S*). (44) 
The gain matrix for this system would be 
G=[1 0.1194 ]. 
The closed-loop eigenvalues could then be computed from 
(A-—BG)V = EVI, 
where 


A, 0 
Vel ve |, r=| 5 ih 
where , is the i” eigenvalue and v; its corresponding eigenvector. In solving the Riccati equation 
the Frobenius norm of the relative residual matrix is 3.14x 10~'°. A simulation result of the stabilized 


system is shown in Figure 12. The initial condition is the Initial Condition Set 1. & 
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Figure 11: Derivative of Lyapunov, ¢ = —0.8 for (a) and (b), {e = 1.7, K = 100} for (c) and (d), 
and {e = 0.08, K = 100} for (e) and (f) 
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279 order singularly perturbed system, stabilized origin 
€= 0.08 


Figure 12: System with stabilized origin (a) phase plane, (b) x,(t), and (c) xo(t) ¢ = 0.08 


Figure 13 (a) and (b) are the plots of the cases of Figure 11 (a)-(b) and (e)-(f) respectively. Initial 
conditions were the Initial Conditions Set 1 for 13 (a), and Initial Conditions Set 3 for Figure 13 (b). 
When € < 0 the system is unstable. Figure 14 (a) and (b) are when kK < 0 in u = —Ksgns. From 
Figure 14 (a) there is a domain of attraction outside of which the trajectory does not converge. 
From simulations with different values of parameters it can be seen that the phase planes when 
€ > 0 while K > 0 in u = —Ksgns have the same characteristic which is represented by Figure 
14 (b). The higher value of the feedback gain K only changes the value of x2 which the trajectory 
converges onto, but the phase plane characteristic is still the same. 


Initial Condition Set 3 The set of n points described by 
n=0 


for i = -50, -30,..., 50 
xi(0) = i, x2(0) =i 
x1(0) = i, x2(0) = -i 


n = n+2 
endfor 
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24 order singularly perturbed system 
€ = -0.8, with variable structure feedback 


2" order singularly perturbed system 


€ = .08, K=100, with variable structure feedback 
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Figure 13: State planes (a) ¢ = 0.3, (b) « = 0.08 
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2" order singularly perturbed system nd : 
‘ ee 2™ order singularly perturbed system 
A S220. B with variable structure feedback (positive 12.) € = .08, K=-100 (positive f.b.), with variable structure feedback 
T T T T T 
r r r r r 
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Figure 14: State planes (a) « = —0.8, (b) ¢ = 0.08 K = —100, ie positive variable structure 
feedback 
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Another second order singularly perturbed system 


Consider the problem 


which can be represented as state equations 


Ly = vr (46) 
EX = —X%1,—-—Lo+u. (47) 


Some knowledge about the system 


The only equilibrium point is the origin. (f(x) = 0 => x, = x2 = 0) The Jacobian matrix is 
OF a OF | 
Se |) cede gy 
x=(«1,0) € € 
1 


~ Oa 
(-1 da/I = 4e). 


The Eigenvalues are shown in Figure 15 graphically as real parts and absolute imaginary part. And 


A 


The Eigenvalues are therefore 


/ 
A 
5- Z | 
= oe 
>] er: 
£ ahs et ee 
no] 
§ 0b 4 
2) 
oa 7 
—5L Z 4 
/ 
/ 
-10 1 | 1 i i | i i i 
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 
Epsilon 
‘| T TS aoa 
0.8 4 
__0.6F A 
£ 
0.47 4 
0.2- a 
0 1 1 1 i | | 1 i 1 
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Figure 15: Real and imaginary parts of the eigenvalues 


the nature of the trajectories can be outlined in Table 2. 
Table 3 shows A, P, 1, A2, V1, and v2 for various «. 
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Values of ¢ | Behaviour of trajectory 
—oo<e<0 saddle point 
e=0 system degenerated 
O0<e<0.25 stable node 
€ > 0.25 stable focus 


Table 2: Qualitative behaviour of the trajectories. 


e <0 
(-) e=-l e=-0.1 € = —0.02 
a 01 0 1 0 1 
ime 10 10 50 50 
fs 05 —05 0.455 —0.5 0.4902 —0.5 
-~05 0 -~0.5 0.45 —0.5 0.49 
a —0.618 —0.9161 —0.9808 
oe 1.618 10.9161 50.9808 
0.8507 —0.7374 —0.7139 
M1 —0.5257 0.6755 0.7002 
0.5257 —0.0912 —0.0196 
0.8507 —0.9958 —0.9998 
e>O0 
(-) e=0.01 e= 01 e=1 
F 0 1 4 0 1 
—100 —100 AO iG pees 
p | | 0505 —0.5 0.555 —0.5 (5 =03 
~0.5 0.505 ~0.5 0.55 65:, <a 
My —1.0102 —T.127 —0.5 + 0.8661 
XS —98.9898 —8 873 —0.5 — 0.8661 
0.7035 0.6637 0.6124 — 0.35361 
"1 O:7107 —0.748 0.7071 i 
—0.0101 —0.112 0.6124 + 0.35361 
ee 0.9999 0.9937 —0.7071 4 


Table 3: 


Solution of the Lyapunov 


equation, the eigenvalues and the eigenvectors 
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Response with various values of ¢ 


When u = 0 the time response of state variables and the state planes are shown in Figure 16 and 
Figure 17 respectively. The result shows that the system is stable when ¢ > 0 and is unstable 
otherwise. According to Table 2, Figure 16 (a)-(c) have an equilibriumpoint as a saddle point, 


(d)-(e) are stable nodes, and (f) is a stable focus. 


A sensitivity analysis of the system was carried out as follows. [Kha96] Letting « = 1 be the nominal 


value of € turns Equation 46 and 47 (u = 0) into 
ty = Ly 
LQ = —X%1— £2. 
And the system with ¢ incorporated can be rewritten as 
1 = £2 
1 1 


Lo = —-X%— -L92. 
& E 


The Jacobian matrix of and of are given by 


OF >. 0 1 
Ox 


and 
Of OT 22) 0 
OX. d¢ | 3+2 


respectively. These Jacobian matrice evaluated at the nominal ¢ yield 


af lek, 
a eae iy Sa 


Ox 
_ 0 
nominal 7 a 25 re , 


and 
of 
OX 


The augmented equation is 


SS PN Gigs Se 8G 
by = (MeN, + [AE], a(t) = 0 
Let 


_ _ | = 0x1 O22 at 
amem| | =[% oz2 | 


nominal 


Attp : //straitstimes.asial.com.sg/lee/leeindex html 


(56) 


( analysis to be continued on page 29.) 
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Figure 16: State variables. (a) ¢ = —1, (b) « = —0.1, (c) e = —0.02, (d) ¢ = 0.01, (e) e =0.1 and 


(f)e=l. 
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Phase plane 
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Figure 17: State planes. (a) ¢ = —1, (b) e = —0.1, (c) e = —0.02, (d) e = 0.01, (e) ¢ = 0.1 and 
(fJe=l. 


Technical report #2, presented to Professor kK. Furuta by Kittisak Tiyapan 29 


With a sign function in feedback 


When 
u = sgn(x; + X2) 


the time response of state variables and the state planes are shown in Figure 19 and Figure 20 
respectively. 


(Continued from page 26, sensitivity analysis.) 


the sensitivity estimation part of Equation 56 is verbosely described as 


and Equation 56 becomes 
ay = wr r1(0) = 10 
oe eee 59) 
t4 = —-%3—-%4+4,+2%_ 24(0) = O 


The result of simulating this system is show in Figure 18. Figure 18 (a) differs from Figure 18 (b) 
in the initial conditions given for 7; and x». Both figures show that x, is slightly more sensitive to 
Ae than 2». In other words, a > ae 


Sensi Aes lots, vanational state variables oe ay Ma, vr ots, vari tation 6) state vandbles 
=0x,/08, X,(0) = 1, X(0) = 1 a -1, 
T T T T 


Figure 18: Sensitivity plot of x3 = Oo and x3 => , (a)x,(0) = x2(0) = 1, (b)x1(0) = —1,42(0) = 1 
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Figure 19: State variables. (a) « = —1, (b) « 
(f)e=1. With a sign function feedback. 
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(f) 


—0.1, (c) e = —0.02, (d) « = 0.01, (e) e = 0.1 and 
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Phase plane, sign FB Phase plane, sign FB 
T T 


Figure 20: State planes. (a) « = —1, (b) « = —0.1, (c) « = —0.02, (d) e = 0.01, (e) ¢ = 0.1 and 
(f) © =1. With a sign function feedback. 
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e | V’s plots | phase plane plots | nature of the equilibrium point 
-1 | Figure 21 | Figure 24 (a), (b) saddle point 
0.1 | Figure 22 | Figure 24 (c), (d) stable node 
0.7 | Figure 23 | Figure 24 (e), (f) stable focus 


Table 4: V of various values of € 


Let the input be 


u = —Ksgn s (60) 
where s is a hyperplane and 
$= ch + ia, CER. (61) 
Then 
&§ = C%,+2£2 
er eee Ksgn (cx, + 2) (62) 
E E E 
Choose the Lyapunov function as 
1 
V= a (63) 
Then 
V = 868 
x Hs Ksgn (cx, + & 
= (cv, +22) [en ee ee 2) (64) 
€ E E 
Let 


c= 0:8; and AK = 15.2; 
the derivative of the Lyapunov function (Equation 64) is plotted the way described by Table 4. 
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00 


dV/dt 


50 


Contour of the derivative of Lyapunov function 
e=-1 


A derivative of the Lyapunov function A derivative of the Lyapunov function 
side, e =-1 3-d,e=-1 


Figure 21: V (a) contour plot, (b) and (c) phase plots ¢ = —1 
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Contour of the derivative of Lyapunov function 
€=0.1 


A derivative of the Lyapunov function A derivative of the Lyapunov function 
side, € = 0.1 3-d,e=0.1 
20 r : teases 
= SSS Ane 
wee 
205 
0 sa a ee 
a 
— 
o4 l 

-20 a ~ : 

-20 
-40 L 3S 40 

3 

60 -| 
-60 iM 

-80 4 i J—— 

5 
a | ~100 = 
-5 < 
-100 - x 
5 -4 -3 -2 -1 0 1 2 3 4 5 
x, x. 


Figure 22: V (a) contour plot, (b) and (c) phase plots ¢ = 0.1 
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Contour of the derivative of Lyapunov function 
€=0.7 


IS \ e/ N 
09 ~20. 7a N71 6 vA 
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Figure 23: V (a) contour plot, (b) and (c) phase plots ¢ = 0.7 
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2" order singularly perturbed system 2" order singularly perturbed system 
€ = -1, without input € = -1, with variable structure control 
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2" order singularly perturbed system 2"4 order singularly perturbed system 
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| 
4b 4b x 4 
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i i 
k 
2b 2 4 
1p 1p @ 4 
OF Of 4 
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x x 
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a} ot | j \ | 
4b 4b | ‘ 4 
5 1 1 1 1 1 1 5 L 1 1 L 1 
-8 -6 4 -2 0 2 4 6 8 -6 4 2 0 2 4 6 


Figure 24: Phase plane of the 2-nd order system (a), (b) ¢ = —1; (c), (d) e =0.1; (e), (f) ¢ = 0.7 
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Fourth order singularly perturbed system 


Consider next the singularly perturbed system 


C= Aux Ajoz Byu, (65) 
es = Agx Aogz Bou. (66) 


Some knowledge about the system 


As this is a linear system the eigenvalues and eigenvectors can be obtained from the matrix A. For 
e=l, 


Pe a | 
= a ee i | 
ee ee ee 
it. di 4h 3B 


The four eigenvalues are 
Ay = 3.7321, A, =1.5+ 0.8667, A3=1.5—0.8662, Aq = 0.2679, 


and the eigenvectors are 


0.445 0.3449 + 0.07761 0.3449 — 0.07767 —0.6388 
» = | 70-2225 |, _ | —0.1343+0.59757 |__| —0.1343- 0.59751]  _ | 0.3194 
0.6894 |’? 0.1343 — 0.59752 |? 3 0.1343 +0.59757 |’ 4 0.6701 
0.5265 —0.1053 — 0.3375 i —0.1053 + 0.33751 —0.2025 


The state equations are converted into a transfer function and the root locus plotted in Figure 25. 


ee gh Ag? re? he 3501 x lo" B _ 
s* — 753 + 16s? -— 15s +3 


Ce tia,DeH1 


Re See HY 


The matrix A for a general case where ¢ # 1 is 


A= 


—_ 

mM 1RmM I rR RO 
—_ 

M [OM | RE RS RS 


mje iro | 
alee ino | 


The movement of the four eigenvalues is portrayed in Figure 26. The values of € used for plotting 
this figure are 
O91 01 180-1778. es O00 oc, 87106: 


In the light of Figure 26 the system is always unstable no matter what value ¢ may take since there 
are always eigenvalues with a positive real-part. With very small value of € (|¢| < 1) the change in 
system characteristic is abrupt. When ¢ = 0 the system degenerates. 
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Real part of eigenvalues 
T 


50 T 


-100- 4 


Figure 25: Root locus of the 4"-order system when € 


Real part of eigenvalues 
T T T 


Imaginary part of eigenvalues 
T T T 


50 


-100- 4 


Real part of eigenvalues, closed—up 


2 T T 


a ae ap yp EE | 
i\ 
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| 
Aaa dal 
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sal rae | 
-1y 4 
-1.54 4 
2 \ i i \ \ f \ 
<2 1.5 -1 -0.5 0 0.5 1 15 2 


Figure 26: Movement of eigenvalues with « (a) real part, (b) imaginary part, (c) real part (closed-up) 
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Response with various values of ¢ 


When ¢ = —1 the state variables and the state plane are shown in Figure 27 and Figure 28 
respectively. Here x3 and x4 are z; and 22 respectively. 


80 10 

70 5 

60 0 

50 -5 
= gy 

40 -10 

30 -15 

20 -20 

10 -25 

0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 
t t 

10 10 

0 5 

-10 0 
fap) st 
x x< 

20 -5 

-30 -10 

-40 -15 

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 


Figure 27: State variable, « = —1. 
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Figure 28: 


State plane, ¢ = —1. 
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When ¢ = —0.1 the state variables and the state plane are shown in Figure 29 and Figure 30 
respectively. 


50 1 
10 

40 
9 
* 30 XY 8 
7 

20 
6 
10 5 

0 0.2 04 06 0.8 0 0.2 04 06 0.8 1 
t t 

10 10 
0 5 

-10 
0 

2-20 x 

5 

-30 
_40 -10 
-50 -15 

0 0.2 04 06 0.8 0 0.2 04 06 0.8 1 


Figure 29: State variable, « = —0.1. 
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x1 
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~20 -20 
-20 -10 0 10 20 -20 -10 0 10 20 


Figure 30: State plane, ¢ = —0.1. 
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When ¢ = —0.01 the state variables and the state plane are shown in Figure 31 and Figure 32 
respectively. 


45 10.2 
40 10.1 
35 10 
30 9.9 
71 X 
25 9.8 
20 9.7 
15 9.6 
10 9.5 
02 04 06 £08 0 02 04 O06 08 1 
t t 
10 10 
0 5 
: | 
0 
2-20 x 


Figure 31: State variable, ¢ = —0.01. 
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Figure 32: State plane, ¢ = —0.01. 
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When ¢ = —0.001 the state variables and the state plane are shown in Figure 33 and Figure 34 
respectively. 
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40 
10.02 
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7; X 
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Figure 33: State variable, « = —0.001. 
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Figure 34: State plane, ¢ = —0.001. 
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When ¢ = 0.01 the state variables and the state plane are shown in Figure 35 and Figure 36 
respectively. 


10 9 


x10 x 10 
4 2 
0 
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2 XQ -4 
-6 
1 
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0 -10 
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t t 
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1 10 
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5 3 
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SQ 2 
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7 1 
1 
0 0 
0 0.02 0.04 0.06 0.08 0.1 0 0.02 0.04 0.06 0.08 0.1 


Figure 35: State variable, ¢ = 0.01. 
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Figure 36: State plane, ¢ = 0.01. 
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When ¢ = 0.1 the state variables and the state plane are shown in Figure 37 and Figure 38 
respectively. 


35 10 
30 9 
25 8 
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0 0.02 004 006 008 0.1 0 0.02 004 006 0.08 0.1 


Figure 37: State variable, ¢ = 0.1. 
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Figure 38: State plane, ¢ = 0.1. 
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When ¢ = 1 the state variables and the state plane are shown in Figure 39 and Figure 40 respectively. 
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400 “ee 
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Figure 39: State variable, « = 1. 
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Figure 40: State plane, ¢ = 1. 
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With a sign function in feedback 


Next use a sign function of the state variables as feedbacks, ie. 
u=sen (41+ 224+ 214+ 22). 


The state planes are shown in Figure 41. The response is unstable. 


Figure 41: State planes when input a sign function of state variables, (a) «= —1, (b)eé=1. 
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In fact the input should have had a negative instead of a positive sign. Equation 65 and 66 turn 
into 


Ly 2° “Q° She - 7 Ly 1 
x —1 1-1 1 x 1 
hes = 9 1 2 1 3 + 1 ‘sgn (xy + 42+ 2+ 22) (67) 
Ce ee em 
22 a aes Fe “Suh 0 


The new system has a structure which varies accordingly with the value of the sign function. Rewrite 
Equation 67 as 
== A=+ Bseny = TE, (68) 


where the switching hyperplane is represented by 
XHM%y4+%24+ 2% 4 2 


and the rest of the symbols are 


v1 2211 1 

_ |2 = eee 1 
a= os ode cee ee ee ee Bl 

te.“ ay 5 A 

22 a 2a 0 


and 
2+sgny 2+sgnyxy Il1+sgnxy Il1+sgny 
—l+sgny 1+sgny —l+sgny 1+sgny 


T= 2+sgn x 1l+sgn x 2+sgn x l+sgn x 
i i i 3 
€ € € € 


The characteristic of the new system depends on eigenvalues of the matrix YT which may take either 
one of these three values, 


Ty xX > 0 
Pad A, y= (69) 
Mey M, 0 
where 

3.3 2 2 1 1 0 0 
0 2 0 2 —2 0 -2 0 
EY Cae Ge 1 ge TS TO GE ae a ng 
a Om iy a> GB 


When xy = 0 the eigenvalues had been shown earlier in Figure 26. The rest two cases are shown in 
Figure 42. Figure 42 (a)-(c) are the case where y > 0 while Figure 42 (d)-(f) are the case where 
x < 0. Figure 42 (d) shows the system to be stable, in other words all the eigenvalues ); satisfies 
the condition 

Re ri < 0, 


if € satisfies the condition 
-l<e<-05, xy <0. 


When xy > 0 the system is always unstable. 
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Real part of eigenvalue, x > 0 Imaginary part of eigenvalues, x > 0 
T T T 
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Real part of eigenvalue, x < 0 
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80 


Figure 42: Movement of eigenvalues with « when x > 0 ((a)-(c)), and when y < 0 ((d)-(f)) 


56 Furuta Lab., Tokyo Institute of Technology 


Variable structure feedback 


Let 
u = —Ksgn s, 


where 


S = CX, + Co + C321 + 2 


is a switching hyperplane. It follows that 


S = CyLy + Coleg + C321 + 22 


The state equations Equation 65 and 66 become 


By = 24, + 294+ 21 + 2g — Ksgn (121 + Coe + C321 + 22) 

L = —%, +249 —2,+ 2 — Ksgn (C121 F Co%q + C321 4 Zo) 

é 2 oD) 2 22 

i + —24 4 — —sen (C121 + Co%2 + 321 + 22) 
€ € € € € 


20 


Let the Lyapunov function be 


The derivative of the Lyapunov function is 


VSS 
1% ZY 2 
= (C121 F Co%q + C32 4 22) + —Zo4 
EE € im 
Cg (—41 + Lo — 21 + 29 — Ksgn (e141 + cote + 321 + 22)) 4 


Cy (241 + 2x + 21 + Zq — Ksgn (cr + C22 + €321 + 22)) 4 


2 ep) 2 22 K 
C3 | 701 Paice Ben (Cutan Cate Caz 2) 
E ee Ss 


1 
= FS [(ciz1 CoX2 + C321 4 22) (1 + 2C3 + 2c1€ = C9€) e+ 


Lo + C3X%_q + 2C1EX2 + C9EXQ + 21 + 263.21 + ClEZ, — CoEZ1 + 


222 + 329 + C1EZq + C2EZ2 — (C3 + (Cy + Co) €) Ken (C121 + C2%2 + 6321 


(78) 


(79) 


The following figures plot derivative of the Lyapunov function with respect to 71, 2, 21, and 29. 


The points plotted follows points shown in the Initial Conditions Set 4. The 


Initial Condition Set 4 k = 0, k is the time step 
for p= -5, -4, ..., 5 
for q = -5, -4, ..., 5 
for r= -5, -4, ..., 5 
for s = -5, -4, ..., 5 
x(k+1). = (p,q; ry °8)5 
endfor; 
endfor; 
endfor; 
endfor; 
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The parameters were 
Cy = Lee C2 = 0.5, C3 = 4.3. 


The feedback gain K is 35 for Figure 43, 44, 45, and 46. These four figures are projections of a four 
dimentional hyper-space onto each of the four axes. 
Figure 43 shows the plot when ¢ = 2. The system is stable since 


V 20VeeDte RR. =e 2D 25.0 CR. 


Derivative of the Lyapunov function Derivative of the Lyapunov function 


e=2 e=2 
0 0 
-500- 4 -500- q 
-1000- aa -1000 4 


-1500- 4 -1500 5 
> > = si 
-2000 -2000 =| 
-2500- 4 -2500 
-3000 =} -3000 
3500 f f ! f ! ! ! f f _3500 f f ! ! 1 f 1 ! f 
- -4 -3 -2 -1 0 1 2 3 4 5 = -4 -3 -2 -1 0 1 2 3 4 5 


(b) 


Derivative of the Lyapunov function 
en2 


Figure 43: Derivative of Lyapunov function (4'"-order singularly perturbed system), ¢ = 2 
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When ¢ = 0.1 the result is shown in Figure 44. A closed-up of Figure 44 shows that the system is 
unstable since dx such that V = 350 > 0. 


Derivative of the Lyapunov function Derivative of the Lyapunov function 
x10° e=0.1 x 10° e=0.1 
0.5 T 0.5 T 


(a) (b) 


Derivative of the Lyapunov function Derivative of the Lyapunov function 


x10° e=01 x 10° e=0.1 
0.5 T 0.5 T 


Figure 44: Derivative of Lyapunov function (4'"-order singularly perturbed system), ¢ = 0.1 


Stability of a system can be determined by using either 
1. the Lyapunov’s stability theorem ({Kha96], Theorem 3.1) or 
2. the Lyapunov’s indirect method ([Kha96], Theorem 3.7.) 


The first one is based on the Lyapunov’s function V () while the second one on whether the Jacobian 
matrix A = OF a) evaluated at the equilibrium point is a stability matrix. Both method has got its 
own unique , namely 


(to be continued on page 59) 
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For ¢ < 0 similar plots can be shown as the plots of Figure 45 and 46 Figure 45 shows that the 
result when ¢ = —1 is unstable. 


Derivative of the Lyapunov function Derivative of the Lyapunov function 
e=-1 e=-1 


1000 


1000 
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| 
1000+ 4 =1000F 


-1500 


(a) 


Derivative of the Lyapunov function 
e=-1 


Figure 45: Derivative of Lyapunov function (4"-order singularly perturbed system), ¢ = —1 


(Continued from page 58, about stability) 


1. the ’s stability theorem only provides a sufficient condition for stability; failure of one Lya- 
punov function candidate does not imply failure of another Lyapunov function (therefore if 
no successful candidate can be found then nothing may be said about the stability or the 
asymptotic stability of the system), 


2. the Lyapunov’s indirect method does not give any about the stability of a system if 


RerA; <0 WV and ai such that Re A; = 0. 


(to be continued on page 60) 
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Figure 46 shows the result when ¢ = —0.01 
a smaller in magnitude of ¢ < 0. 


Derivative of the Lyapunov function 
x10 e=-0.01 
3 T 
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. The result shows that the system is also unstable for 


Derivative of the Lyapunov function 
e=-0.01 


25 
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0.5 L L L L Ll Ll L L Ll 
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Derivative of the Lyapunov function 
e=-0.01 
3 T 


ob 
ae 


(b) 


Derivative of the Lyapunov function 
e=-0.01 


3 T 


Figure 


(Continued from page 59, about stability) 
The Lyapunov 


46: Derivative of Lyapunov function (4"-order singularly perturbed system), ¢ = —0.01 


PA+A™P=-Q, Qpositive definite (80) 


is the link between the two method. In order to find out whether A is a stability matrix or not, one 
can either 


1. find all the eigenvalues A; of A, if Re A; < 0 then A is a stability matrix, or 


2. solve Equation 80, if it has a definite solution P then A is a stability matrix. 
(to be continued on page 61) 
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Next increase the value of the gain K. In comparison to Figure 44 increase the gain K ten times 
(ie K is now 350) and show the result in Figure 47. From Figure 47 the system became stable after 
the gain was increased. 


Derivative of the Lyapunov function Derivative of the Lyapunov function 


x 10° e=0.1,K=350 x10° e=0.1,K=350 
0 0 
1b 4 1k 4 
2b 4 2b 4 
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4b 4 S 
\ \ ZA 
Ne ENG 
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5b 4 
6 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 
-5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5 
x x 
1 2 
(a) (b) 
Derivative of the Lyapunov function Derivative of the Lyapunov function 
5 €=0.1,K = 350 5 £=0.1,K =350 


Figure 47: Derivative of Lyapunov function (4"-order singularly perturbed system), ¢ = 0.1, K was 
increased ten times. 


(Continued from page 60, about stability) 
Which method is better depends on whether there is any uncertainty in the system . That is, 
1. if A is a constant matrix, then finding the eigenvalues is better recommended since ;’s also 


give direct information about the response of the linearized system. 


2. if A is a perturbed matrix, then finding the solution of the Lyapunov equation is necessary 
because its solution leads to a Lyapunov function, the existence of which allows conclusions 
about the stability of the system. » 
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However in the case when ¢ = —1 < 0 increasing the gain Kk ten times did not make the stable 
stable, compare Figure 48 with Figure 45. 


Derivative of the Lyapunov function Derivative of the Lyapunov function 


25 T T T T T T T T T 2.5 


0.57 4 0.5- 


(b) 


Derivative of the Lyapunov function 
e=-1,K=350 


Figure 48: Derivative of Lyapunov function (4'"-order singularly perturbed system), ¢ = —1, K was 
increased ten times. 


When a Lyapunov function V(x) satisfies the conditions 


V(0) =0, V(x) >0Vx Ee D—O, and (V)(x) <0 VeEeD, DCR, 
x = 0 being the equilibrium point, it does not mean that D is the domain of attraction for the 
system because there is no guarantee that a trajectory which starts from x(0) € D will remain 
forever in D. It only means that the equilibrium point is stable and, if furthermore 


V(r) <0, Vere D—-O 


then the equilibrium point is asymptotically stable. The domain on which the Lyapunov function 
is valid has got nothing to do with the domain of attraction of the system. ~) 
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Phase plane plots of the system 


A couple of phase plane plots are shown for 
1. ¢ =2, K = 35 (Figure 49 (a) and (b)) which is stable, 


2. ¢=0.1, K = 35 (Figure 49 (c) and (d)) which is unstable except perhaps within a limitted 
domain of attraction, 


3. € =0.1, K = 350 (Figure 49 (e) and (f)) which is stable. 


The parameters were 
Cy = soe C2 = 0.5, C3 = 4.3. 


And the initial conditions used were the Initial Condition Set 5. 
Initial Condition Set 5 The following set points 


e five points on the line connecting (—7, —6, —4, —13) and (8,5, 12,3), namely 
(-7,-6, —4, -13), (—4, -3.8, 0.8, -9.8), (—1, 1.6, 2.4, -6.6), (2,0.6, 5.6, 3.4), and 
(5, 2.8, 8.8, 0.2) 


e five points on the line connecting (14, 7,5, —5) and (—3,5,—10,4), namely 
(—3, -7, 10, —5), (0.4, —4.6, —7, 3.2), (3.8, -2.2, 4, -1.4), (7.2,0.2,—1,0.4), and 
(10.6, 2.6, 2, 2.2) 


e five points on the line connecting (—6, 13,4, —3) and (21, —3, —4,8), namely 
(—6, —3, 4, -3), (0.6, 0.2, -2.4, —0.8), (4.8, 3.4, 0.8, 1.4), (10.2, 6.6, 0.8, 3.6), and 
(15.6, 9.8, 2.4, 5.8). 


e five points on the line connecting (—4,6,5,4) and (5, —3, —5, —2), namely 
(—4, -3, -5, 2), (—2.2, -1.2, -3, -0.8), (—0.4, 0.6, -1, 0.4), (1.4, 2.4, 1, 1.6), and 
(3.2, 4.2, 3, 2.8). 


e five points on the line connecting (—3, —8,3,—1) and (14, 2,—1,5), namely 
(=8, 8,1, —1); (0.4, =6, —0.2,.0.2), (3:8, —4.0,0.6, 1,4), (7:2; -2, 1:4, 2.6), and 
(10.6, 0, 2.2, 3.8). 


e five points on the line connecting (—10,7,—8,—12) and (13, —9.3, 7, 1.2), namely 
(—10, —9.3, —8, —12), (—5.4, —6.04, —5, —9.36), (—0.8, —2.78, —2, —6.72), 
(3.8, 0.48, 1, —4.08), and (8.4, 3.74, 4, -1.44). 


Figure 49 (a)-(f) show these plots. Comparing Figure 49 (c)-(d) with Figure 49 (a)-(b), when 
varepsilon becomes smaller the system becomes less stable. Note that the plot of the Lyapunov’s 
derivative (Figure 43 and Figure 44) also shows that ¢ = 0.1 is less stable than ¢ = 2 since Figure 43 
is negative-semidefinite while Figure 44 is neither negative-definite nor negative-semidefinite. The 
stable region of Figure 49 (a) and (b) (¢ = 2) is larger than that of Figure 49 (c)-(d) (¢ = 0.1.) 
When the feedback gain AK was increased (Figure 49 (e)-(f)) the system becomes more stable. The 
simulation time was t < 3sec. 
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Figure 49: Phase planes (4'"-order singularly perturbed system) (a)-(b) © = 2,K 


e=0.1,K =35, (e)-(f) ¢ =0.1,K = 350 


20 T 
157 4 
107 4 
5b 
8 
BS) 
~ oF Boren SSS 4 
o SPR, 
5h oe oo 4 
oC 2 
~10- % 4 
15+ 4 
_20 1 1 1 1 1 1 
-20 -15 -10 -5 0 5 10 15 20 
Zz 
1 
4'"_order singularly perturbed system 
Zo2 en OAK = 38) : 
20 T T 
157 a 
107 =| 
5} | | 
oO 
~ oF KK 4 
ox & 
2 eae, 
5b 4 
@ 
~10F 4 
15+ 4 
_20 1 1 1 1 1 1 
-20 -15 -10 -5 0 5 10 15 20 
Zz 
1 
4"_order singularly perturbed system 
2,-Z, (€ = 0.1, K = 350) 
10 
————~9 
5L | 
x—_____»—0 
=—=,°? 
x—9 %59 
ob —& & Gace 4 
OO gtx 
° es } 
na er 
o—__ S$ x—9 
x—O 
5b oe # 
Q—_x 8 
10h oy 4 
-15 1 1 1 1 1 
-10 -5 0 5 10 


35, (c)-(d) 


Technical report #2, presented to Professor K. Furuta by Kittisak Tiyapan 65 


A polygonal phase plane 


This section is an attempt to find a way of displaying movement of more than three state variables 
on a two dimensional plane. This method is named a polygonal phase plane has the following 
description. 


1. The axes for the state variables are represented by lines going outward from the origin. The 
axes are equally spaced around the origin. 


2. Each axes has a normallized scale from r = 0.1 for the minimum value of the simulated state 
variable to r = 1 for the maximum value, r is the radius starting from the origin. 


3. Plot the state variables as snapshots of polygons each having a vertex on each axis. The 
number of snapshots is < 10. Also plot the starting point as a circle and the end point as a 
triangle. 


4. The order of the axes starts from the one which goes toward the north as the first one or 271 
with the the other state variables having an axis in a direction on a line rotated clockwise. 


5. Whether the system converges or not can be seen from the simulated numerical data. If the 
data is still varying with each step at end of the simulation time, then it is probable that it 
diverges. 


Figure 50 (a) to (f) show polygonal phase plots of the system. Following is the description of each 
subfigure. 


Figure 50 (a) is when the parameters were 
€ =0.01, c, = 3.7, co = 1.5, cz = 0.3, K = 500. 


The initial conditions was (7, —8, 4,3). The range of the axes are normalized with the original 
values of [6.742, 7| for x1, [—8.2553, —8] for x2, |—22.2527, 4] for z,, and [2.7901, 3.0648] for zo. 
The simulation time was less than 3 seconds. The system is stable and the result ends on a 
hyperplane. 


Figure 50 (b) is when the parameters were 
E= 0.1, C= 0.7, cg = 1.5, C3. = 2.3, K = 50. 


The initial conditions was (3, —4,6,0). The range of the axes are normalized with the original 
values of [2.8527, 3] for 71, [—4.1956, —4] for x2, [4.6955, 6] for z1, and [0, 0.1358] for zo. The 
simulation time was less than 3 seconds. The system is stable and the result ends on a 
hyperplane. 


Figure 50 (c) is when the parameters were 
e=1, cy = 2.4, co = 1.2, c3 = 2, K=5. 


The initial conditions was (5,1,3,2). The range of the axes are normalized with the original 
values of [5,1.1 x 10°] for x1, [—49950, 0] for x2, [3, 1.51 x 10°] for 21, and [2,1.25 x 10°] for 
zg. The simulation time was less than 3 seconds. The system is not stable and the result did 
not converge. 
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Figure 50 (d) is when the parameters were 
620 ee Te SOS = a Ra. 


The initial conditions was (—7,—2,—6,1). The range of the axes are normalized with the 
original values of [—7,6.37 x 10°] for x1, [—2, 4.63 x 103] for x2, [—9.89 x 10°, —6] for 21, and 
[—510, 7.4] for z2. The simulation time was 3 seconds. The result did not converge therefore 
the system is unstable. 


Figure 50 (e) is when the parameters were 
PS GeNH LT. es 1S 8. K = A008. 


The initial conditions was (0,4, 12,7). The range of the axes are normalized with the original 
values of |—4.03, 0] for 71, [—0.199, 4] for xo, [8.08, 12] for 21, and [7, 7.26] for zo. The simulation 
time was less than 3 seconds. The system converges. This can be compared with Figure 50 
(c) which also has ¢ = 1 but which has a smaller value of K. 


Figure 50 (f) is when the parameters were 
e=0.001, 4 = 19,6 = 1.3, c= 258, K-= 400, 


The initial conditions was (—5,0.3,7,—3). The range of the axes are normalized with the 
original values of [—5.0014, —5] for x1, [0.2986, 0.3] for x2, [5.5883, 7] for z1, and [—3.0156, —3] 
for z2. The simulated time was less than 3 seconds. Because there are very few snapshots in 
the result, one knows that the system converged quickly onto a hyperplane. 


The following idea is new, and has yet got to be formally realized. I will do this in the course of 
the making the subsequent technical reports. The brief outline is drawn as follows. 


1. A general polygonal phase plane for high dimension system may be drawn (in similar manner 
to the results given above) by increasing the number of axis. Generally speaking, if a system 
is n-dimensional with x;,7 = 1,2,...,n as its state variables, then the angles these state 
variables make clock-wise with respect to the vertical 0-o’clock axis are 


Ok 
pee. P2019. cae 
nr 


2. A normalized polygonal phase plane may be drawn with all the axis normalized to the nominal 
value of the state variable it represents. This idea of normalized values is used extensively in 
power system engineering. One possible drawback could be the fact that, unlike in the field of 
power system engineering where the parameters does not vary dramatically away from their 
nominal values, the state variables in control system engineering could fluctuate with virtually 
no limit. 


3. With the aid of computer softwares the polygonal phase plane plots could be put together to 
form a movie showing a continuous movement of the trajectory. 
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4. In order to distinguish the new trajectory in an n-dimensional polygonal phase space from 
the trajectory of a conventional state plane, a new terminology adopted from other existing 
disciplines of science metamorphosis shall be used to describe the movement of the phase 
polygon shown by these snapshots. This terminology (ie phase-metamorphosis) will be better 
equipped to accommodate a video sequence (a it phase-metamorphosis video) made from these 
snapshots (a phase-metamorphosis diagram.) 


5. Lastly, (conventional) phase planes (2-dimensional,) phase spaces (3-dimensional) as well 
as the time-domain plots of the state variables may be incorporated into polygonal phase- 
metamorphosis diagram as small separate boxes. 


This outline will be improved upon in other technical reports. 


ore) 
The degenerate system 
Then let ¢ = 0 and obtain a degenerate system 
29 = Ajo®+ Bou, 20) = 2(0) (81) 
22 = Apx®+ Bou, 2°(0) #4 (0) 
where 
Aj = Au _ Ay2A55 Ao1 (82) 
BY i By = Aj A55 Bo (83) 
Ay. Se Ae Asi (84) 
By = —Agy By (85) 
The matrix Ag: is assumed to be nonsingular. The effect of degeneration is two folds, ie 
1. reduce the state variable z from its status as state variables 
2. desert the initial condition z(0) of z. 
; 1-1-3333 
From the values of Aj, Ajo, By, Aoi, Ao2 and By above it follows that A, = 0 1 : 


o _ | 0.6667 _ | -1 —0.3333 
om 2 a= 0 —0, 3333 
The system reduces in dimension from 4“"-order to 
two eigenvalues are both 1’s. 

A descriptive diagram comparing the structure of the original equations (Equation 65 and 66) with 
that of the degenerated system (Equation 81) can be seen in reference ({Nai88], Figure 2.3.) 


| ana 2 =| 


gnd 


—0.6667 
0.3333 
-order. The new system is unstable since the 
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Polygonal phase plane Folygonal phase plane 
e=0.01, C,=3.7, c,=1.5, c,=0.3, K=500 e=0.1, c= 0.7, C= 1.5, ¢. g=e-3, K=50 


<> 


(a) (b) 


Polygonal phase plane ‘enue aie plane 
oy .2, C,=2, K=! e=-0.1, c,=2.1, K=32 


(c) (d) 


goicake phase plane Polygonal phase Hane 
e=1,¢,=1 c,=2.3, K=400 e=0.001, c, =1.9, ,=1.3, Cy =2.8, K=400 


(e) (f) 


Figure 50: Polygonal phase planes (a) ¢ = 0.01, (b) ¢ = 0.1, (c)e = 1, (d)e = -0.1, (e)e = 1 
(with a higher gain), and (f) « = 0.001 
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The state variables and the state plane are shown in Figure 51 and Figure 52 respectively. (Errata: 
The labels x1, 22, x3 and x4 written on the plots are x9, x8, z? and 29 respectively.) 
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Figure 51: State variable, degenerated system. 
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Figure 52: State plane, degenerated system. 
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The degenerate system with a sign function feedback 


Input a sign function of state variables to the degenerated system and obtain Figure 53 (a). Add 
state feedbacks as well as the sign function input and obtain Figure 53 (b). (these simulations may 
be incorrect due to error in programming, I think it might have been done on a 4“"-order system 


by mistake) 
50 50 \ 


x nO 
-50 -50. \ 
-50 50 -50 0 50 
x1 x2 
50 50 
_—S—_ 
4 °F O———[—[— = 
SSS 
-50 -50. 
-50 0 50 -50 0 50 
z1 x1 
(a) 
50 50 


x NO 
-50 -50. 
50 50 -50 0 50 
x1 x2 
50 50 


Figure 53: State plane, degenerated system with (a) u = sgn (a, +%24+2%+ 2%), (b) u = 
sgn (41 + %2 + 2% + 22) 
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Degenerated system with variable structure control 


In this section the origin of the degenerated system was 


1. stabilized via state feedback having a gain obtained from the solution of the Riccati equation 
described in page 19. This amounts to Equation 86 with K = 0. The result is shown in Figure 
55 (a), then 


2. a variable structure was then applied by using K 4 0 in Equation 86. Figure 55 (b)-(f) show 
the response. In Figure 55 (f) Equation 87 with 7 = 10 was used instead of Equation 86. 


When simulating with a negative feedback of sign function (K < 0) the numerical iteration suffers 
indefinite repetition once the trajectory reaches the surface of the hyperplane. Therefore with 
negative feedback Equation 87 was used as an approximate of Equation 86 in simulating. An 
interesting observation when simulating with kK < 0,y = 10 in Equation 87 is that when c > 0.99 
the trajectory converges toward the origin (eg. see Figure 54 (a),) but when c < 0.98 the system 
converges to two equilibrium points which are away from the origin (eg. see Figure 54 (b).) 

The control input can be written in a general form as 


u = —2.8229 x, — 1.8227. + K sen (cx, +X), (86) 
and a special control used as the approximate for Equation 86 
u = —2.8229 x, — 1.822%. + K tanh y (cr, + 22). (87) 
Solving Equation 42 (a Riccati equation) for this system leads to 


p= 3.4844 0.2499 _ | —1.43804 
~ | 0.2499 0.8281 ~ | —2.0972 


| , G= | 2.8229 1.8228 |, 


where [ is a matrix containing the closed-loop eigenvalues and G is the gain matrix. The Frobenius 
norm of the relative residual matrix is 6.836 x 107!°.. 


Degenerated si Higa patio syoltn wth. VS co control rated si ingularly pert oad eye atv c8 control 
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Figure 54: Phase plane (a) K = —8,c=1.5 (b) K = —-8,c=0.3 
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Degenerated singularly perturbed system 
the remaining two state—variables 


Degenerated singularly perturbed system with VS control 
30,c=1.7 


15+ 
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(c) 
Degenerated singularly perturbed system with VS control 
K =30,c =-1.7 
20 T 


Degenerated singularly perturbed system with VS control 
=3,c=17 


(d) 


Degenerated singularly perturbed system with VS control 
with tanh instead of sign fn., K = -30, c = 1.7 


Figure 55? Phase plane-(a) ik =0,(b) Kh = 376 = 7, (eo) i = B06 = LT, (ed) Sj eS — LF, 


(e) K = 30,c = -1.7, and (f) K = —30,c=1.7 
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Boundary value problem 


29" April 1998 
Consider a boundary value problem [Nai88] 


d?x dz 
= 88 
“ae | a oe 
with boundary conditions 
(Oost (89) 
The solution of this equation is 
g295 bbs 
CS + 4 90 
ee —1 és my 


was shown in Figure 56. Figure 56 (d) shows the process of degeneration (ie. changing from ¢« 4 0 
to ¢ = 0) which is opposite to the process of singular perturbation (ie. changing from ¢ = 0 to 
é # 0) and which is characterized by the loss of some boundary conditions. Figure 56 (g) is the 
combination of Figure 56 (a) to (f). 


Another boundary value problem 


Next consider another boundary value problem |[Nai88] 


dex’ 4 dz 
ea tel we —(0 91 
dt? e dt % 2h) 
with boundary conditions 
HO) = 15... ei a3 (92) 


The solution of this system is 


a [-13 ERR AMES 4 3 4 1g (oe )'-13) 


el/2 (eeev eae G. vi-te = a) Ga wists + 1) -1 ea (13 el/2 vinfe_1/2e71 = 3) 


ot/2 CERE 41/2e71 1/2 UV) Guam - 1) Ga vi=te sp i) (03) 


Figure 57 shows the time response of this system. 
Figure 57 (d) shows the case where degeneration occurred, resulting in the loss of one boundary 
condition. Figure 57 (g) is Figure 57 from (a) to (f) represented in one graph. 
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Figure 56: Boundary Value problem, (a) « = —0.1, (b) e = —0.05, (c) « = —0.03, (d) « = 0.03, (e) 
€ = 0.05, (f) « = 0.1, and (g) combined in a single plot 
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Figure 57: Boundary Value problem, (a) ¢ = —0.1, (b) e = —0.05, (c) e = —0.03, (d) e = 0.03, (e) 
é = 0.05, (f) ¢ =0.1, and (g) combined into a single plot 
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E (—00, -}] (-4, 1) [1, co) 
equilibrium point is a node focus node 


Table 5: Dependence of characteristic of the trajectory upon € 


Perturbed system 


24" April 1998 
Consider the system 


gy) = —%, +a(1+a1)+e(1+21)’ (94) 
—21(41 +1) +4. (95) 


£9 


To find the equilibrium point, let f(x) = 0 (where « = f(z,u), and u = 0) and solve the simulta- 
neous equations 


—2, +49(1+2,)+e(1+2)? = 0 (96) 
—21(r, + 1) (97) 


The equilibrium point is at (0, —e). 
The Jacobian matrix is 


Osi. —142e(l+a1)+%. l+a 
Ox SIS 2a Ole - 
This Jacobian matrix when evaluated at the equilibrium point gives 
Of —-1l+3.¢ 1 | 
pee = (99) 
OX |,-(0,-c) 7 ie 
Eigenvalues of A are 
1 3 1 
Mg==5 Gets Vo bet Oe, (100) 


The qualitative behaviour of the trajectory is summerized in Table 5 depends on the value of ¢, for 
example 


—3-66+9e7 >0 (101) 


the trajectory has a node which may either be stable or unstable depends on whether A; < 0 or 
> 0 respectively. Moreover when 


a 
3 
the trajectory becomes structurally unstable since the eigenvalues on the imaginary axis is vulner- 
able to perturbations. The equilibrium point becomes a hyperbolic equilibrium point. 
Figure 58 shows plots versus ¢ of real and imaginary parts of the eigenvalues. 
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Eigenvalues 
Real and Imaginary parts 
4 T T 


Real part 


Imaginary part 


Figure 58: Eigenvalues, real and imaginary parts 


Domain of attraction 


The domain of attraction can be found whenever the matrix A in Equation 99 is a stability matrix 
(ie all of its eigenvalues are positive definite or Re; < 0.) The plot of Re); of the system in Equation 
94 and 95 when u = 0 is shown in Figure 59 (a) and (b). 

As an example, consider a case when ¢ = 0.1, which means 


SS) Sesto) — S07 i 
PE ae 9 


Solving the Lyapunov equation having Q = J, 


PA+A?TP=-Q=-I (102) 
and get as a result 
a ce oe | ee iss | 
Then the Lyapunov function can be chosen as 
V = (1.4286 2; + 0.522) 21 + (0.521 + 1.7786 x2) 22. (103) 


And its derivative is 


Vi = —0.71428 03 + 2? (—3.8572 — 0.6 22) + 
v1 (—1.51089 + x2) (—0.189107 + 22) + #9 (0.14 22). (104) 
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Real part of the eigenvalues 
Real part of the eigenvalues close-up 
T 


Re 
Red 
° 
] 


-\ 


Figure 59: Real part of eigenvalues (a), and close-up (b) 


0.28572 2, — 3.8572 27 — 0.71428 x? + 0.1 a2 — 1.72122 —0.6 2222+ 43+2,23 (105) 
—3.8572||2||2 + 0.29||z||2 + 0.72||x||3 + 0.1||z||2 + V1.72? + 1||x||2 + 


i 
(Fle) (vO Tl) (106) 
—1.4949]|a||2 + 1.8862||x|/3 <0, for ||a|/3 =r? < 0.7925. (107) 


IA 


Thus V (a) is negative definite in a ball ([Kha96],Chapter 4) 


D = {x € R? | |I2\2 <r} (108) 


of radius 


r = V0.7925 = 0.8903. 


The simplest estimate is the set 


2, = {2 € R?|V(a) <ch. (109) 
Then 

ae er Vig) SA Pir (110) 
leads to 


c=0.8 <_ 1.0738 x 0.7925 = 0.8510. 


Then better estimates of the region of attraction can be obtained by finding the largest ball on 
which V(x) is negative definite. This can be done via the trajectory reversing method ([Kha96], 
Chapter 4.2.) This method is essentially done by backward integrating the state equations using 
points which lie on the surface of the boundary of the estimate of the region of attraction. Backward 
integration of the state equations can be done by forward integrating the new state equations 


é =—f(2). (111) 
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Let wu = 0 and write Equation 94 and 95 in the form of Equation 111 


Ly = Vy — £2(14+ 21) =e(1+2)? (112) 
ai (ay +1): (113) 


Lp 


Then simulate these new equations using the 20 obtained by solving 


20k 
0, = ao" and (114) 
T 
9 | cosOz cos 6;, - 
Pk sin 0; | < sin 0; | =e (115) 


as initial conditions. 

The contours in Figure 60 are the estimated region of attraction obtained by the method of 
trajectory-reversing described above. Backward integration was done over the interval [0,t;], 7 = 
1,2,3,...,44. At t; = 0.4seconds the simulation encountered singularities and the further simula- 
tion is no longer possible. The wedge-like shape was due to backward integration of three points 
on Wo which approach infinity rapidly. 


Estimates of the region of attraction 
trajectory-reversing method, k = 44 
1 T T T T T i 


0.5 


0.55 


Figure 60: Estimates of the region of attraction (¢ = —1) 
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Sensitivity analysis of the system 


81 


The procedure used to do the stability analysis is referred to a reference [Kha96]. Let the nominal 
value of perturbing parameter be ¢ = 1. The nominal system is then given by (u = 0) 


ry 


Lg 


—%4 ro(1 L1) (1 


= —21(r1 + 1). 


The Jacobian matrix of was given by Equation 98 and OL is 


These Jacobian matrice when evaluated at the nominal value of € give 


of 
Ox 


nominal 


ala Ol aa) sey Lewy Of 
, and 
—l1— 2%, 0 


Of _ (Pen)? 
=| 0 | 


Oe 


Augmented equations are ([Kha96], Equation 2.11) 


Let 


Ti 
LQ 
£3 
La 


GS Fle AY: x (to) 
Poa as [286490] | Eee as) 


—%,4 x2(1 L1) (1 oa ae 


—21 (21 + 1), 


nominal | 


1 
2 


(—1+2(1+41) +22) 23+ (1+21) 24+ (1+2)*, x3 


(-1 = 221) X3, 


Figure 61: Sensitivity plot 
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Response with various values of ¢ with no input 

Phase planes of this system (u = 0) are shown in Figure 62 (¢ = —1.5), Figure 63 (« = —1), Figure 
E 


64 (« = —0.1), Figure 65 (« = 0), Figure 66 (¢ = 0.3), Figure 67 (¢ = 1.3) and Figure 68 (¢ = 1.7). 


Phase plane 
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Figure 62: Phase plane. ¢ = —1.5 


Technical report #2, presented to Professor K. Furuta by Kittisak Tiyapan 83 


x2 


x2 


Phase plane 


Figure 63: Phase plane. ¢ = —1 
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Figure 64: Phase plane. ¢ = —0.1 
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Phase plane 
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Figure 65: Phase plane. ¢ = 0 
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Figure 66: Phase plane. ¢ = 0.3 
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Phase plane 
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Figure 67: Phase plane. ¢ = 1.3 
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Phase plane 
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Figure 68: Phase plane. ¢ = 1.7 


With a sign function in feedback 

Introducing u = sgn(x, + x2) as an input and obtain Figure 69 (¢ = —1.5), Figure 70 (« = —1), 
Figure 71 (« = —0.1), Figure 72 (« = 0), Figure 73 (€ = 0.3), Figure 74 (e« = 1.3) and Figure 75 
(C= 157) 
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Figure 69: Phase plane. ¢ = —1.5 with VS input. 
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Phase plane 
10 T 


x2 


Figure 70: Phase plane. ¢ = —1 with VS input. 
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Phase plane 
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Figure 71: Phase plane. ¢ = —0.1 with VS input. 
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Figure 72: Phase plane. ¢ =0 with VS input. 
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Phase plane 
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Figure 73: Phase plane. ¢ = 0.3 with VS input. 
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Figure 74: Phase plane. ¢ = 1.3 with VS input. 
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Phase plane 
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| | | | 
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Figure 75: Phase plane. ¢ = 1.7 with VS input. 
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Variable structure control 


Let the input be 
u = —Ksgn (cv, + 22), (123) 


where K,c € R and the switching hyperplane is s = cx, + 29. 
And let the Lyapunov function be 


1 
= 5° 

Then 

$= —(1+(1+21)) +¢(-x e(1 7) + (14 x1) £2) — Ksgn (cx, + £2) (124) 
and 

V = s8 

= (cy +2) (—(1+ (L+m)) +¢(—a1 te (1+?) +(1 +21) 2) - 
Ksgn (cx, + £2)) (125) 


One could plot Equation 125 for various values of ¢ (as has been done in dmy earlier works,) in 
order to find out whether this system is stable. But at this moment let us look at another possiblity 
of analysing a variable structure control. 


Variable structure control with a hyperbolic tangent function 


This section introduces and investigates how a hyperbolic tangent function could be used to estimate 
the sign function. One great advantage of the hyperbolic tangent function over the sign function is 
that it can be differentiated, as 
dtanh(z) 
dx 


Figure 76 shows a juxtaposition of the tanh(yx) and sgn (x), y = 1,5, 10,50, and 100. The figure 
shows that the hyperbolic tangent function can be used to approximate the sign function very well. 


= sech(z). 


Let the input to the system is 
u = K tanh (y (cv, + 22)). (126) 


Equation 126 is an approximation of Equation 123. 
Equation 94 and 95 with Equation 126 as an input would turn into 


gq = -2-4 r2(1 r1) e(1 x1)" (127) 
fg = —21(4, +1) — K tanh (y(cr1 + 22)), (128) 


where 7y,c, K are constants. 
The Jacobian matrix is then 


le. GLO): 4 —142¢(1+21)+ 22 1 ey (129) 
~ Og | —1 — 2x1 — cyKsech? (y (cx; + t2)) —yKsech? (y (cr, + 22)) |” 


If the location of the equilbrium point is known and one would like to find a domain of attraction, 
one could use the following procedure 
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tanh(x) and sign(x) 
solid line is tanh(x), dashed line is sign(x) 
t 


Figure 76: Hyperbolic tangent function compared with sign function 
1. evaluate the Jacobian matrix at the equilibrium point to find find the matrix A of the system 
linearized around the equilibrium point, 
2. solve the Lyapunov equation for P, 


3. obtain the Lyapunov function from 
V =e Px, (130) 


4. find the first estimate of the domain of attraction by looking at terms in the Lyapunov function, 
and 


5. plot better estimates of the region of attraction by the trajectory-reversing method. 
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But in our case instead of solving the Equation 127 and 128 which could prove to be quite difficult, 
I will simplify the problem by letting 
cal. “kh 30s “=50,- and! .c=-1,3; 


and plot the response in state space as Figure 77. From Figure 77 the response is stable. 


Perturbed system with variable structure control 
with tanh instead of sign fn., K = 30, y= 50, € = 0.1,c=1.3 


5 0 5 19 15 20 25 


Figure 77: Hyperbolic tangent variable structure control, state plane 


Since in Figure 77 the equilibrium point is in the vicinity of the origin it is possible to locate it by 
calculating f;(x) and f2(x) and look for a point where these two functions are at minimum value. 
The calculation done using Algorithm 1 shows that the equilibrium point is at the origin and 


I ealateans = Ara r2(1 x1) e(1 wa)? =e: 


| 131 
Fl enialibrinn == —21 (a1 F 1) — K tanh (y (cary + X2)) = 0. ( ) 


Algorithm 1 fi_eq = 10; f£2_eq = 10; x1_eq = 0, x2_eq = 0; 
for x1 = -5 upto 5 step 0.1 
for x2 = -5 upto 5 step 0.1 
calculate f1(x) and f2(x); 
if (lf1(x)| < fi_eq) and (|f2(x)| < £2_eq) 
xl_eq = x1; x2_éq = x2; flleq = f1(x); f226eq = £2(x); 
endif 
endfor 
endfor 


Do Algorithm 1 again for various values of ¢ then plot Figure 78, which shows that 


(11, %2)|equitibrium = (0,0) for ¢<0.27 (kK = 30,7 =50,and c = 1.3). 
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xX, and x, 
when min(f, (x) and 1,0) 


f, and f, 
when min(f,(x) and £304)) 
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Figure 78: (a), (c), x1 and £2; (b), (d) fi(x) and fo(x) when min (fi(x)) and min (fo(z)) 


The Jacobian matrix at the equilibrium point for ¢ = 0.1 < 0.27 is then 


-—1 1 
a —1951 —1500 | 
and the eigenvalues are correspondingly 
Mt = —2.3, and v2 = —1498.7. 


Then solve the Lyapunov equation AP + PA? = —I to get P which is 


p—| 90-2175 —0.2825 
~ | —0.2825 0.3677 |’ 


and is a positive definite matrix since the eigenvalues are all positive (0.003 and 0.5849.) 
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Let the Lyapunov function be 
V=2'Pr= (0.2175 x, — 0.2825 xg) x1 + (—0.2825 x1 + 0.3677 x2) x2 
Then the derivative of Lyapunov function is 
Vo = (0.2175 21 — 0.2825 x9) #1; + x; (0.2175 & — 0.2825 9) + £p (—0.2825 4, + 0.3677 &) + 
(—0.2825 x1 + 0.3677 x2) & (132) 
= (0.21752; — 0.2825 x2) (—a, + 0.1(1 + 21)” + (#1) a2) + 


x, (0.2175 (—x1 + 0.1 (1 + a1)” + (1 + a) @2) — 0.2825 (— (#1 (1 + 21) - 


30 tanh (50 (1.3.1 + w2)))) + a2 (—0.2825 (—a1 + 0.1 (1 + a1)”) + (1+ 21) 22) + 

0.3677 (— (a, (1+ 71)) — 30 tanh (50 (1.32, + %2)))) + 

(—0.2825 x1 + 0.3677 x2) (— (21 (1 + 21)) — 30 tanh (50 (1.321 + 22))) (133) 
= 0.6085 x? + x? (0.217 — 0.3569 v2) + (—0.0565 — 0.565 x2) v2 — 

0.565 x, (—0.442364 + x2) (0.174045 + x2) + (16.95 7 — 

22.062 x2) tanh (65 x; + 5022) (134) 
= 0.0435 2, + 0.217 x? + 0.6085 22 — 0.0565 x2 + 0.1516 2122 — 0.3569 x72. — 

0.565 23 — 0.565 2123 + 16.95 x, tanh (50 (1.3 2, + r2)) — 

22.062 x2 tanh (50 (1.3 21 + x2)) (135) 


Equation 134 is in a simplified form and Equation 135 in an expanded form. 
In polar coordinates 


x, = psiné, and (136) 
tg = pcosé. (137) 


The Lyapunov function in Equation was derived from state equations which are written in x; — %9 
Cartesian coordinates. Change of the variables x, and x2 in Equation into a p—@ polar-coordinates 
results in 


V(p,0) = p cos(6) (0.3677 p cos(@) — 0.2825p sin(@)) + p sin(@) (—0.2825 p cos(@)+ 
0.2175 p sin(@)) (138) 
p’ (0.3677 cos?(0) — 0.565 cos(6) sin(9) + 0.2175 sin?(6)) . (139) 
Similarly, the derivative of the Lyapunov function in Equation 132 to 135 becomes 
V(p,0) (—0.2825 p cos(@) + 0.2175 p sin(6)) (—p sin(@) + p cos(6) (1 + p sin(6)) + 
0.1 (1+ p sin(9))”) + p sin(6) (0.2175 (—p sin() + p cos(9) (1 + p sin()) + 
0.1(1 + p sin(6)) *) — 0.2825 (—p sin(@) (1 + p sin(@)) — 30 tanh (50 (p cos(@)+ 
1.3 p sin(9))))) + p cos(@) (—0.2825 (—p sin(@) + p cos(@) (1 + p sin(@)) + 
( 
) 


0.1(1 + p sin(6)) *) + 0.3677 (—p sin(@ (1+ p sin(@)) — 30 tanh (50 (p cos(@)+ 

1.3 p sin(@))))) + (0.3677 p cos(@) — 0.2825 p sin(@)) (—p sin(@) (1+ p sin(@)) — 

30 tanh (50 (p cos(@) + 1.3 psin(@)))) (140) 
=p a (1 + p sin(@)) + cos(0) (—0.0565 + 0.1516 p sin(#)— 

0.3569 p sin?()) + sin(@) (0.0435 + 0.217 p sin() + 0.6085 p? sin?(0)) + 

(—22.062 cos(@) + 16.95 sin(@)) tanh (50 p cos(@) + 65 p sin(A))) (141) 
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Notice that both the Lyapunov function and its derivative was derived in the Cartesian coordinates 
%1,%_2 domain. The change of variable was only adopted to facilitate the consideration of the 
equations. This is different from deriving a new Lyapunov function from the polar-coordinates p, 


domain whose states are 
r = \a2+23 (142) 


x 
6 = tan}. (143) 
pal 
Following is shown a wrong method in deriving the state equations in polar-coordinate, then the 
right method is given afterward. 


The wrong method is to differentiate Equation 142, 143 and then substitute x, « in the right hand 
side of the equation obtained. To illustrate that this method does not work, try solving for state 
equations with parameters as before. 


d 
eae xt + x 
224 X4 + 2% 9X9 
ee (144) 
2\/xt + x3 
: d oD) 
é = —tan>— 
dt XY 
228 of 29 
= a (145) 
14+3 
Substitute Equation 136, 137 as well as 
é, = sin0p+cosé pb (146) 
é = cosOp—psindd (147) 
into Equation 144 and 145 to obtain 
p 2p sin(6)(sin(0)@ + cos(@)p0) + 2.cos(@)p(cos(0)p — p sin(0)@) 
2\/cos?(6) p? + p?sin?(0) 
=; 
biz = cot (0) csc(6) (sin(0) A + cos(@)p0) + csc(0)(cos(0)~ — p sin(6)9) 
7 p(1 + cot?(6)) 
22, 
Which is obviously wrong! tit 


The correct method is as follows. 
From Equation 127 and 128 substitute Equation 136, 137, 146, and 147 to obtain two simultaneous 
equations 


sin(@) p+ cos(4) 90 —(psin(0)) + cos(@)p(1 + p sin(@)) +e (1+ p sin(@))? (148) 


cos(@)p — p sin(0)@9 = —psin(@)(1 +p sin(@)) — K tanh(y(cos(@)p + cp sin(@))) (149) 
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which are then solved for p and 0. 
The state equations are then 


p ~ cos%(6) F sin?(6) (-e sin(#) + p sin?(0) — 2ep sin?(0) — ep? sin®(0)+ 
K cos(@) tanh (y(p cos(@) + cp sin(6)))) (150) 
= sin(@)(e+ (—1+ 2€)p sin(@) + ep sin(@)) — K cos(@) tanh(yp(cos(6) + ¢ sin(@))) (151) 
: 1 2 
Cos = poet Oye paar) (cos(@ ) (0 sin(@) — p cos(@) (1 +p sin(@) —e(1+ p sin(6)) ) — 
sin(@) (p sin(@)(1 + p sin(@)) + & tanh (y(p cos(@) + cp sin(@)))))) (152) 
= a (c (4 + p) cos(@) + p (4 — ep cos(3@) + 4p sin(@) — 2 sin(20) + 4e sin(20)) + 
4K sin(@) tanh (yp (cos(@) + c sin(@)))) . (153) 


fit 


Now let us look at Equation 141 and try to find the first estimate of the region of attraction. From 
Equation 141 


V 


= —0.565 p* cos*(@) — 0.565 p® (cos? (0) sin(0)) + 
p cos(@) (—0.0565 + 0.1516 p sin(0) — 0.3569 p” sin?(0)) a 
p sin(9) (0.0435 + 0.217 p sin() + 0.6085 p? sin?(0)) + 
p (—22.062 cos() + 16.95 sin(@)) tanh (50 p cos(@) + 65 p sin(@)) . (154) 


By noticing that 


—0.565 p” cos?(#) <0, and tanh (50: cos(@) + 65 pp sin(@)) < 1, 


Equation 154 becomes 


V< 


IA 


= 


= 


p (—0.0565 cos(@) + 0.0435 sin(@) — 22.062 cos(@) + 16.95 sin(@)) 

-p” (—0.565 + 0.1516 cos(@) sin() + 0.217 sin?(6)) 

+p (—0.565 cos?(0) sin(@) — 0.3569 cos(@) sin?(0) + 0.6085 sin*(0)) (155) 
p (| — 22.1185] - |cos(@)| + |16.9935] - | sin(@)|) 

+p (—0.565 + |0.1516| - | cos() sin(6)| + [0.217 - |sin?(0)}) 

+p* (| — 0.565] - |cos?(8) sin(0)| + | — 0.3569 - |cos(0) sin?()] + [0.6085] - |sin*(6)|)(156) 
p (22.1185 x 1+ 16.9935 x 1) + p? (—0.565 + 0.1516 x 0.5 + 0.217 x 1) 

+p? (0.565 x 0.3849 + 0.3569 x 0.3849 + 0.6085 x 1) (157) 
31.112 p — 0.267 p? + 0.9633 p®> << 0. (158) 


Equation 158 does not yield any solution p € R therefore the estimate of the region of attraction 
has not yet been obtained. 
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Trying to find the estimate from Equation 135 only leads to another failure since from Equation 
135 


; 0.1516 
V < (0.0435? + 0.05652 + V16.952 + 22.0622) ||x||2 + (0.2172 + 0.5652 + Ilr /3 + 
2 


V0.3569? + 0.5652\ . 13 
llr llo 


; (159) 


< 0.9428]|x|/3 + 0.681||a||3 + 27.8928]|2||2 (160) 


(0.008 


which also yields no solution ||z||p € FR. 

This could be because our assumption that the origin is the equilibrium point was wrong. Looking 
at Equation 131 at the origin f(z) = ¢, instead of zero. The Lyapunov function used above 
was obtained from solving the Lyapunov equation using the matrix A obtained with the above 
assumption. So if the assumption was not sound, A loses its validity. 

Let us look more carefully at the response of the system. Plot Figure 79 with the Initial Condition 
Set 6 which is similar to the Initial Condition Set 1 often refered to before but with a smaller scale. 
Figure 79 shows that the system has a stable focus at approximately (0.5, —1). 


Initial Condition Set 6 The set of n points described by 


n=0 
for 1 = -0.5 to 0.5 step 0.02 
n = nti, calculate the n_th initial values as x1i(0) = i, x2(0) = 
n = n+1, calculate the n_th initial values as x1(0) i, x2(0) 
endfor 


| 
H- 


i} 
| 
H- 


Let us do the simulation again to investigate the behaviour of the system inside the limit cycle by 
shifting the Initial Condition Set 6 downward along the x-axis and obtain the Initial Condision Set 
7. Use the new initial condition set to simulate Figure 80. From the response shown the trajectory 
starting from inside the limit cycle moves outward toward it. 


Initial Condition Set 7 The set of n points described by 


n= 0 

for i = -0.5 to 0.5 step 0.02 
n = nti, calculate the n_th initial values as x1(0) =i, x2(0) = i-1; 
n = n+1, calculate the n_th initial values as x1(0) =i, x2(0) = -i-1; 


endfor 
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Perturbed system with variable structure control, state plane 
with tanh instead of sign fn., K = 30, y= 50, €=0.1,c =1.3 
T T 
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Perturbed system with variable structure control, x,(0) Perturbed system with variable structure control, x,(t) 


with tanh instead of sign fn., K = 30, y= 50,¢=0.1,c=1.3 wath tanh instead of sign Ins = 30; y=S0,e = O:1;¢54.9 
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Figure 79: Tanh variable structure controlled response (a) state plane, (b) x(t), (c) x2(t) 
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Perturbed system with variable structure control, state plane 
with tanh instead of sign fn., K = 30, y= 50, € = 0.1,c=1.3 
T T T ‘f T T 
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Figure 80: Tanh variable structure controlled response (a) state plane, (b) x,(t), (c) x2(t) inside 
the limit cycle 
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With the light of Figure 79 and Figure 80 one could use, say a radius 0.5 circle centered at (0.5, —1), 
as an ad hoc first estimate of the region of attraction boundary, and then proceed on with the 
trajectory-reversing method to find better estimates. 
But let us try an approach which is a little bit more systematic. Let us unperturbed the system 
by setting « = 0, then use the matrix A obtained to find the region of attraction. When ¢ = 0 
Equation 127 and 128 yield the origin as an equilibrium point. The Jacobian matrix becomes 

Of (x) oe Nee Let 


= 161 
Ox —1— 22%, —cyK sech (y (c41 + 42)) —yK sech (y (cx + 22)) len) 


and at the equilibrium point becomes 


-1 1 
Me ap (162) 


With the values of parameters used above 


-1 1 
oa —1951 —1500 | a8) 


The eigenvalues of this matrix are 


Solve the Lyapunov equation (Equation 102) for P using A in Equation 163 leads to 


0.2175 —0.2825 
pe —0.2825 0.3677 | eo 
hence the Lyapunov function can be found from Equation 130 to be 
V = 0.2175 x? — 0.565 212 + 0.3677 x3 (165) 


which the derivative is 


V 0.13 2? + 0.565 x? + 0.2646 x12 — 0.3004 222 — 0.565 22 — 0.565 2123 + 


16.95 x; tanh (50 (1.3.71 + x2)) — 22.062 x2 tanh (50 (1.3 2; + x2)) (166) 


< (0.13 2} + 0.2646 x22 — 0.565 3) + (0.565 x} — 0.3004 xpr2— 
0.565 r123) + (16.95 x; — 22.062 x2) tanh (50 (1.3 x; + r2)) 
< 16.95? + 22.062?||2||> + (0.13 + 0.2646 — 0.565)||a||3 + 


(167) 


0.30042 + 0.5652 
(0.50 vi ag ) iets 


2 


Equation 167 does not yield a solution ||z||2 € R either. 

Since we could not find a first estimate of the region of attraction, then the trajectory-reversing 
method can not be applied. But I will try using the trajectory-reversing method by starting from 
Figure 79 and 80 which show a limit cycle close to 


(x, — 0.5)? + (a2 +1)? = 0.57. (168) 
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Starting from 20 points on the contour of Equation 168 using the convention given in page ‘points’ 
simulate the inverse system in the form given by Equation 111. The Equation 127 and 128 with 
the set of parameter values used above become 


ty — #2 (1+2,) —0.1(14+ 24) (169) 
x,(1+2,)+ K tanh (7 (cr, + 22)). (170) 


4 


i) 
And the plot upto the 49” time step, ie 
0 <t; < 0.01 seconds, Ve 2 si A 
is shown in Figure 81 (a) while in Figure 81 (b) the time steps were 
0 <t,; < 0.5 seconds, P= 195 Or. 


Figure 81 (c) and (d) are both closed-ups for Figure 81 (b). 

The figure shows a trajectory which seems to be on a fixed course, which could be a characteristic of 
a system subjected to a control input. Eventhough my approach in obtaining Figure 81 comes from 
looking at a phase plane response (therefore ad hoc) instead of at the Lyapunov function (which 
would have been more mathematical, ie more intrigue) I do believe that this approach works with 
a system if it is continuous. 


The following is only to fill in the blank space of this and the next pages. 
Three inequalities often used in considering the upperbound of a Lyapunov function are 


lvi] < lala, (171) 

1 
leit) < 5|lellz, and (172) 
lax, tbrq| < Var+b?\IzIlo, (173) 


where ||||2 = \/27 + x3. 


About a hyperbolic tangent function, written in exponential form it is 


tanhxz = See (174) 
Ca e* 


Its taylor series upto the 10°’-order is 


ge. De ATE! .. bog? 
tanha = 2 — — 4 + Ol (x). 175 
ara ae (we yO (2) ne) 


And the taylor series of 
f = tanh (y(cr1 + %2)), 
where y,c € R and upto the 5“’-order is 
58 
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sec, (b) t < 0.5 sec, (c) and (d) closed-up’s of (b) 
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A 4 order singularly perturbed system 


This section is the addition and correction to the result in Technical Report #1 |Tiy98a] (see also 
[Tiy98b]). Consider the same system to a system studied in Technical Report #1, 


g = Arv+Bu (177) 
= Czxr+Du (178) 
with 
—8 —)5 1 1.25 
4 0 0 0 
ca 0 QO -l11 -—2.5 |’ 
0 O 4 0 
0 
0 
B=! 9 99843 |? 
0 
Co— —2.82843 1.59099 0.70711 0.88388 
and 


D=0. 


Response of the system without an input 


Let e = 1 


Incorporation of ¢ 


The system equations becomes 


¢ = EAr+E'Bu (179) 
= Cr+ Du. (180) 
when 
E500 Q0 
Oy ol AQ" 0 
a OS as 1 1) 
00 06 


is incorporated into the system equations. Notice that here E~'B = B and E~'A only differs from 
A in the 4" row which has only got one element (column 3.) Also CE~'B = 2 regardless of the 
values of €. 


Some knowledge about the new system 


The equilibrium point still remains at the origin but the eigenvalues are affected by the ¢ introduced 


and are eo 
—44+)2, and — oe aa Bee 
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Phase plane plots of the new system with different ¢’s 


The following figures are the plots of the system without a control input. When ¢ = 1 the state 
plane plot were shown in Figure 82. The initial conditions were taken as 


Initial Condition Set 8 7s a set of points described by 
for i=-5,5 


for j=-5,-3,-1,1,3,5 
for (p,q)=(1,2),(1,3),(1,4), (2,3), (2,4), (3,4) 


x_p(0O) =i 
x_q(0O) = j 
x_k(0O) = 0, for all k which are not equal to p or q 
endfor 
endfor 
endfor 


Figure 83 shows the state planes obtained when ¢ = 0.01. 
Figure 84 shows the state planes obtained when ¢ = 10%. 


State variables plot 


The following are the time domain plots of the state variables for the values of €’s previously used. 
Figure 85 shows the state variables and state metamorphosis when ¢ = 1. The phase metamorphosis 
diagram was simulated upto t = 6 seconds. Shown on the diagram, the term ftgna) is used to refer 
to the time when all the state variables converge to the origin to within +0.001 tolerance. Notice 
that tanal < 6 seconds The initial conditions used were 


Initial Condition Set 9 is a set of point described by 


for i= -3, 3 
for j = -3, 3 
for k = -3, 3 
for m = -3, 3 
x_1(0) = i, x_2(0) = j, x_3(0) =k, x_4(0) =m 
endfor 
endfor 
endfor 
endfor 


Figure 86 shows the state variables and state metamorphosis when ¢ = 0.01. 
Figure 87 shows the state variables and the metamorphosis of states when ¢ = 107. 
Figure 88 shows the result simulated for ¢ = 0.1 using Initial Condition Set 9 as the intial conditions. 
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A singularly perturbed system, without control 
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A singularly perturbed system, without control A singularly perturbed system, without control 
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Figure 83: State diagrams of a system with « = 0.01 
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Figure 84: State diagrams of a system with « = 10~+ 


109 


110 
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state variable x,(t), eal 
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Singularly perturbed system, no input 
state variable x(t), e=l 
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Figure 85: (a) to (d) state variables of a system with « = 1, (e) phase metamorphosis diagram 
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Singularly perturbed system, no input Singularly perturbed system, no input 
state variable xX, th e€=0.01 state variable Xft €=0.01 
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Figure 86: (a) to (d) state variables of a system with « = 0.01, (e) polygonal phase metamorphosis 
diagram 
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Figure 87: (a) to (d) state variables of a system with « = 10~*, (e) polygonal phase metamorphosis 
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Singularly perturbed system, without input 
output y(t) 
20 T T T T 


Singularly perturbed system, without input 
state plane Xs-Xy 


Figure 88: Output y(t) (a) and state diagrams ((b) and (c)) of the system with « = 0.1, no input 
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Response of the system with a VS feedback input 


Design a controller with an input signal 


u = Cz +sgn(CE7'Ax + CE'Bu), (181) 
where 
thi, a et) 
son(f) =< 0). - f= 0 
—-1 ,f <0. 


Equation 181 could not be simulated because algebraic loops were encountered (according to MAT- 
LAB’s error messages.) Therefore modified equation used instead was 


u = Ca +sgn(CE'Ax). (182) 


When ¢ = 1 the result can be shown as in Figure 89. Figure 89 is quite similar to Figure 82 where 
the system was without any control and with the same ¢ = 1. Here the initial conditions used were 
the same as those used for Figure 82. 


The next picture, Figure 90, shows the result of a simulation when ¢ = 0.01. Figure 90 should be 
compared with Figure 83 which is the case with the same value of ¢ (¢ = 0.01), but which has no 
control input. Both Figure 90 and Figure 83 has the same set of initial conditions. 


For the case where ¢ = 10~ the result is show in Figure 91. Figure 91 can be compared with Figure 
84 to see the effect which the feedback applied has on the system response. The initial conditions 
used in Figure 91 were the same as those used for obtaining Figure 84. 


The plots of the state variables when a variable structure feedback is applied to the system with 
€ = 1 are shown in Figure 92 (together with a phase metamorphosis.) The phase metamorphosis 
diagram shown in Figure 92 (e) has (1,1,1,1) as the initial condition. 


When ¢ = 0.01 the state variables plot in time domain are shown in Figure 93. 


The state variables plot in time domain when ¢ = 107+ are shown in Figure 94. In Figure 92, 
93, and 94 a hyperbolic tangent function was used instead of a sign function, so that instead of 
Equation 181 the input is 


u = Cx + tanh (70(CE-'Ar + CE“'Bu)). (183) 
Then for comparison with Figure 88, Figure 95 was plotted. 


Discussion 


The two poles which are not affected by the value of ¢ (the ones located at —4 + 21.) have the 
natural frequency w,, = 4.4721 and the damping factor 0.8944. (The natural frequency of a linear 
time invariant system pole is 

__ [1og() 


i = 
where t, is the sampling time. The damping factor is 


¢ = —cos(/(log ))), 


n 


A being the eigenvalue. 
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= Se ee 


4 


(a) 


A singularly perturbed system, with VS feedback 
e=1 


-6 L 


5 -4 


A singularly perturbed system, with VS feedback 
e=1 


Furuta Lab., Tokyo Institute of Technology 


A singularly perturbed system, with VS feedback 
e=1 
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A singularly perturbed system, with VS feedback 
e=1 
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Figure 89: State diagrams of a feedback system with « = 1 
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Figure 90: State diagrams of a feedback system with ¢ = 0.01 


Lae 


118 


AN 


8.0 MAA tp tt 
$§ 4 3 2 4 #O 4 2 8 4 5 
x, 


Furuta Lab., Tokyo Institute of Technology 


Figure 91: State diagrams of a feedback system with « = 10~* 
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Figure 92: (a) and (b) state variables of the system with V.S. feedback, ¢ = 1, (e) metamorphosis 


when controlled with Tanh V.S. 
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Singularly perturbed system, with VS input Singularly perturbed system, with ys input 
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Figure 93: (a) and (b) state variables of the system with V.S. feedback, ¢ = 0.01, (e) metamorphosis 
when controlled with Tanh V.S. 


Technical report #2, presented to Professor kK. Furuta by Kittisak Tiyapan 121 


Singularly perturbed system, with ys input Singularly perturbed system, with ys input 
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Figure 94: (a) and (b) state variables of the system with V.S. feedback, ¢ = 10~*, (e) metamorphosis 
when controlled with Tanh V.S. 
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Singularly perturbed system, with VS feedback 
output y(t) against t 
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Figure 95: Output y(t) (a) and state diagrams ((b) and (c)) of the feedback system with « = 0.1 


